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Abstract 



We study nonequilibrium dynamics of SU(2) pure gauge theory 
starting from initial over-population, where intense classical gauge 
fields are characterized by a single momentum scale Q s . Classical- 
statistical lattice simulations indicate a quick evolution towards an ap- 
j^H' proximate scaling behavior with exponent 3/2 at intermediate times. 

ch , Remarkably, the value for the scaling exponent may be understood 

as arising from the leading 0(g 2 ) contribution in the presence of a 
time-dependent background field. The phenomenon is associated to 
weak wave turbulence describing an energy cascade towards higher 
\^0 , momenta. This particular aspect is very similar to what is observed 

t^J" ' for scalar theories, where an effective cubic interaction arises because 

VO , of the presence of a time-dependent Bose condensate. 
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O ■ 1 Motivation and overview 



The question of thermalization is one of the most pressing issues in the 
physics of ultra-relativistic heavy- ion collisions. These experiments involve 
S^ | far-from-equilibrium systems of strongly interacting matter described by 

quantum chromodynamics (QCD). The data reveal robust collective phe- 
nomena, whose theoretical understanding from QCD represents a great chal- 
lenge pp. This concerns even the theoretically "clean case" of very large 
nuclei and energies. In this idealized limit, the matter formed shortly after 
the collision may be described by classical gluon fields with a characteristic 
"saturation" momentum scale Q s that grows with the energy and the size 
of nuclei. Since Q s is large, the strong coupling constant a(Q s ) is small. 
However, the gluons at saturation are strongly correlated because their low 
momentum modes have very high occupancy ~ 1/a. 

Recently, it has been argued that in the earlier stages of a heavy-ion 
collision at sufficiently high energies, the parametrically large gluon density 



may lead to the phenomenon of Bose condensation [2j. This has to be com- 
pared to alternative scenarios [3J. For a recent review see Ref. [3]. Inelastic, 
particle number changing processes preclude the possibility that the true 
equilibrium state be a Bose condensate, but there is the possibility that a 
transient condensate develops during the evolution of the system. Finally, 
this question is related to the competing time scales for elastic as compared 
to inelastic processes in gauge theories. However, the weak-coupling anal- 
ysis is complicated by the fact that modes with occupancies ~ 1/a require 
non-perturbative descriptions. 

Important aspects of this non-perturbative gluon dynamics can be de- 
scribed by classical-statistical lattice gauge theory simulation techniques 
El [71 E]. Such classical field simulations, with suitable averages over 
the initial conditions, can accurately describe the dynamics if the expecta- 
tion values of field anti-commutators are much larger than the corresponding 
commutators [9]. Loosely speaking, this is realized in the presence of suffi- 
ciently high occupation numbers per mode, n(p). This concerns regimes of 
large occupation (n(p) 3> 1) all the way up to the over-populated regime, 
where n{p) ~ 1/a. They fail of course when typical occupations become 
smaller than unity, at which point quantum effects need to be taken into 
account. The range of validity of classical-statistical simulations has been 
tested to high accuracy for self-interacting scalar field theories [101 [HI [12] 
and theories with fermions [13], where appropriate far-from-equilibrium ap- 
proximation techniques are available also directly for the respective quantum 
field theory. 

In particular, the dynamical formation of a Bose condensate was recently 
demonstrated for scalar quantum fields starting from initial over-population 
|14j . It has been shown that Bose condensation occurs as a consequence 
of an inverse particle cascade with a universal power-law spectrum. This 
particle transport towards low momenta is part of a dual cascade, in which 
energy is also transfered by weak wave turbulence towards higher momenta. 
Exponents associated with these two cascades are under analytical control 
and are well reproduced in the simulations [121I15J. In particular, the value 
of the exponent for the UV cascade can be understood as arising from the 
existence of the condensate leading to an effective cubic interaction [16]. 
Condensation was also discussed in Ref. [17] . and similar dynamics was an- 
alyzed in Ref. [18] in the context of cold atoms emphasizing also the role 
of non-trivial topological configurations [19]. Whereas a Bose condensate 
develops and remains present in the relativistic system until eventually in- 
elastic scattering depletes it, for the non-relativistic theory with conserved 
particle number no decay of the condensate due to number changing pro- 



cesses is observed [T] 

In this work, we study the non-abelian dynamics of over-populated gauge 
fields for the SU(2) gauge group. A major focus is to work out further the 
relevant differences and similarities between gauge and scalar degrees of 
freedom out of equilibrium in this context. The most remarkable result will 
be that classical simulations starting from initial over-population indicate 
for Yang-Mills theory the same turbulence exponent as for scalars. This 
appears non-trivial for the following reasons. In scalar theories, the fact 
that the interactions are nearly local in momentum space typically plays 
a crucial role to explain the flow of momenta. This is apparently not the 
case in gauge theories, where one may go from hard to soft momenta in one 
collision |20j. 

Furthermore, for scalar theories the inelastic processes are of higher or- 
der in the coupling than elastic collisions. Therefore, one expects that the 
role of number changing processes, which can compete with condensate for- 
mation, can be suppressed for weak coupling in scalar theories. In contrast, 
elastic and inelastic processes are parametrically of the same order in the 
gauge theory. It is, therefore, a question beyond simple parametric estimates 
whether condensation occurs from initial over-population or not. 

We explain that both scalar and gauge theories can show the same tur- 
bulent scaling exponents in the presence of a time-dependent background 
field. This can be done using resummed perturbation theory, since the 
phenomenon of weak wave turbulence occurs in the kinetic regime with oc- 
cupancies in the range 1/a ^> n(p) ^$> 1 [21] . In this regime one may safely 
choose to study occupancies which, in our case, are derived from equal-time 
correlation functions in Coulomb gauge. 

For the gauge theory it is a subtle question of how to interpret the most 
infrared modes beyond the kinetic regime in this way. It has been argued 
that the infrared cascading solutions observed for scalars [12] can be carried 
over to non-abelian gauge theories [22]. It would be very desirable to devise 
suitable gauge-invariant measures, in particular, of condensation. This non- 
trivial task is beyond the scope of the present work. 

This Letter is organized as follows. In section [2] we describe the classical- 
statistical lattice gauge theory approach, and the results obtained from 
initial over-population. In section [3] we compare the numerical results to 
perturbative computations of scaling exponents in the presence of a time- 
dependent background field. We conclude in section HI 



2 Classical lattice gauge theory simulations 

We study the real-time evolution of classical-statistical Yang-Mills theory 
following closely Ref. [6], to which we refer for further technical details. 
For the simulations we employ the Wilsonian lattice action for SU(2) gauge 
theory in Minkowski space-time: 



S[U] 



-AEEJ^I ( Tr U *S* + Tc Ul fil ) - 1 

x i *■ 

+aEE(^t ( Tr u ** + Tr U U) - l 



X i<j 



\2Trl 



with x = (x°,x) and spatial Lorentz indices i,j = 1,2,3. It is given 



in terms of the plaquette variable U x n U = U x uU x+ n V Jj ,~ u Ux V , where 



Ux,vn = U x ^ v . Here U x „ is the parallel transporter associated with the 
link from the neighboring lattice point x + p, to the point x in the direc- 
tion of the lattice axis \i = 0,1,2,3. The definitions /3q = 2^/Trl/gQ and 
(3 S = 2Trl/(gg7) contain the lattice parameter 7 = a s /dt, where a s de- 
notes the spatial and at the temporal lattice spacings, and we will consider 

9o = 9s = 9- 

Varying the action ([1]) w.r.t. the spatial link variables U x j yields the 
classical lattice equations of motion. Variation w.r.t. to a temporal link 
gives the Gauss constraint. We define the gauge fields as 

gAUx) = -^-Tr(a a U t (x)) (2) 

2a s 

where a 1 , a 1 and a 3 denote the three Pauli matrices. The coupling g can 
be scaled out of the classical equations of motion and we will set g = 1 for 
the simulations. The initial time derivatives A^x = 0, x) are set to zero, 
which implements the Gauss constraint at all times. Shown results are from 
computations on spatial lattices with iV 3 = 128 3 — 256 3 sites. 

In order to make contact with discussions in the literature, which are typ- 
ically formulated in terms of Boltzmann-type equations, we need to extract 
suitably gauge- fixed distribution functions. These are related to equal-time 
correlation functions, which are obtained by repeated numerical integration 
of the classical lattice equations of motion and Monte Carlo sampling of ini- 
tial conditions [6]. The dynamics is solved in temporal axial gauge. Though 
this gauge is very efficient for computational purposes, it involves already 
on the perturbative level spurious poles in the propagator which makes the 



extraction of sensible distribution functions or dispersion relations difficult. 
Therefore, we choose to transform the configurations to Coulomb gauge. The 
gauge transformation to achieve Coulomb gauge is calculated by the overre- 
laxation method described in Ref. [23]. We typically use 10 5 overrelaxation 
steps, which achieves a very good convergence for the employed lattice sizes. 
The Coulomb gauge distribution functions can safely be extracted for mo- 
menta where perturbation theory is useful, although it is maybe unclear 
how to interpret the most infrared modes in this way. Fortunately, we will 
see that important lessons can already be learned for not too low momenta 
such that a characterization of physics in terms of distribution functions is 
useful. 

Distributions will be computed with the help of time-dependent field mo- 
mentum modes, which are obtained as Af(t,p) = f d 3 xAf(t,~x.) exp(ip -x). 
We average over color and Lorentz indices for better statistics. The occu- 
pation number distribution we then define as 

f 



«*(*) = v v(m,p)\ 2 u\A(t,p)\*) cl , (3) 

where (. . .) c i also denotes the classical average over typically five to ten runs 
and V is the volume. Correspondingly, we consider the "dispersion" 



(!£(*, P )l 2 )ci 

which in Coulomb gauge should only deviate from a free, linear behavior in 
the regime of large occupancies at low momenta. 

As mentioned above, we want to study the time-evolution starting from 
initial over-population. This means that at early times the distribution has 
a non-perturbatively large occupancy n p (Q s ) ~ 1/g 2 at the characteristic 
momentum scale Q s . We use Gaussian initial conditions such that all the 
relevant information about the initial state is contained in the (equal-time) 
two-point correlation functions. In Fig. Q] we show the occupation number 
distribution in Coulomb gauge as a function of spatial momentum for differ- 
ent times in units of Q s . The (red) dashed-dotted line shows the approximate 
initial distribution as it is shortly after the electric field components built 
up, since we start from Ef = —Af = at t = to fulfill the Gauss con- 
straint. Relatively quickly we observe the evolution towards a distribution, 
which for an intermediate time range can be well described by a power-law 
for momenta |p| < 2Q S . The shown (blue) dashed curve at tQ s = 315 and 
the (black) solid one at tQ s = 730 indicate that the power-law behavior is 
rather well described by the dashed-dotted fit curve ~ |p|~ 3 ' 2 - 



Q_ 



CO 



0.1 



0.01 



0.001 



0.0001 



tQ s «0 

tQ s =315 

t Q s = 730 

■s. f . tQ s =1575 

\tCL = 4410 
\ b 




0.1 



p/Q £ 



Figure 1: Occupation number distribution as a function of momentum for 
different times. 



Once this power-law behavior is established, the subsequent evolution 
becomes rather slow or quasi-stationary. As time proceeds, at some point 
deviations from a simple power become visible. The (grey) dotted curve for 
tQ s = 1575 shows already a somewhat steeper behavior at lower momenta 
and a diminished slope at higher momenta. At tQ s = 4410 the curve is 
clearly not described by a simple power. In order to characterize the transi- 
tion to and the evolution away from the ~ |p|~ 3 ' 2 behavior, we fit a power- 
law dependence. The resulting exponent as a function of time is shown in 
Fig. [TJ The error bars correspond to the change of the power as the fitting 
region is varied by choosing the lower bound from (0.3 — 0.6)Q S and the up- 
per bound from (0.9 — 1.3)Q S . The results indicate that initially the system 
evolves rapidly to 3/2 (long-dashed curve) and spends a relatively long time 
around that value, before it slowly evolves to lower values. Of course, this 
evolution towards lower values is expected since the system will thermal- 
ize classically to a distribution ~ |p| _1 at sufficiently late times (which in 
general happens rather quickly when the occupied high-momentum modes 
reach the highest available momenta on the lattice). However, the observa- 
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Figure 2: The exponent of a fitted power law to the occupation number dis- 
tribution displayed in Fig. [T]as a function of time. The error bars correspond 
to the fluctuation of the power as the fitting region is varied as explained in 
the main text. 



tion that the intermediate evolution is very slow and stays also a significant 
time around the value 4/3 is consistent with earlier investigations starting 
from initial conditions featuring plasma instabilities [21] , This adds to our 
observation that the observed apparent scaling behavior happens for a wide 
range of initial conditions. 

Before we analyze the power-law behavior in more detail in section [3] 
below, we present in Fig. [3] the dispersion. For sufficiently high momenta 
the considered quantity is always close to the free expression oj p = |p| as it 
should. We find that discrepancies from the massless dispersion are present 
only below p 2 < 0.1Q 2 , with decreasing amplitude as the time grows. Above 
we explained that the occupation number distribution seems to follows a 
power-law for intermediate momenta at not too late times. Fig{T]shows that 
the UV part of the spectrum is very slowly filled up, moving the breakdown 
of the power-law dependence in the direction of high momenta. In the IR, 
the breakdown of the power-law dependence is roughly at the same momenta 
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Figure 3: Dispersion as a function of momentum for different times. 



where the dispersion starts to deviate from the free dispersion. 



3 Perturbative scaling exponents 

For quantum fields A" (x) one can define two independent connected two- 
point correlation functions out of equilibrium, which may be associated to 
the anti-commutator and the commutator, 



F$(X,V) 



1 



({A-(x),Al(y)})-A^x)A b M, 



p<Z(x,y) = i([Al(x),Al(y)}), 



(5) 



respectively. Here we took into account a possible expectation value or 
"background field" 

Al{x) = (A-(x)). (6) 

Loosely speaking, the spectral function p determines which states are avail- 
able, while the statistical propagator F contains the information about how 
often a state is occupied. The spectral function is related to the retarded 




Figure 4: Gluon part of the one-loop contribution to the self-energy with 
(2PI) resummed propagator lines. The crossed circles indicate an effective 
three-vertex in the presence of a background gauge field potential. 

propagator Gi R \ and the advanced one G(a) as p = Gr R \ — Gim . A tremen- 
dous simplification of thermal equilibrium is that the spectral and statistical 
functions are related by the fluctuation-dissipation relation, which is not as- 
sumed here [9]. 

For instance, the one-loop self-energy correction to two-point correlation 
functions in the presence of a background gauge potential is diagrammat- 
ically given in Fig. [H To avoid problems of secularity in standard per- 
turbation theory, here the lines are meant in the two-particle irreducible 
(2PI) effective action scheme, where self-energies are expressed in terms of 
self-consistently dressed propagators [9] . This includes also the background- 
field dependence of the propagators. The crossed circles indicate an effective 
three-gluon vertex gV^, which is appearing at one loop in the presence of 
the background gauge field potential © [24J. This effective three- vertex 
consists of the conventional tree-level vertex and an ^.-dependent term, 

yabc yabc + yabc /^ 

The standard tree-level part reads in Fourier space 

rabc 



Vn 



(p,q,k) = f abc [g^{p-q) 1 J rgu 1 {q-k)^+g 1 ^{k-p) v J , (8) 

where f abc are the structure constants of the non-abelian gauge group. Fi- 
nally, we will be interested in a situation where the background field has 
a residual (space-) time dependence. Then the corresponding part of the 
interaction (J7|) reads in configuration space 

YA%ry( x >Vi z ) = ( C ac,bd9nwA*(x) + C ab)dc g uy A*(x) + C ahfid # 7M A*(x)) 

gS d+1 (x-y)6 d+1 (x-z) (9) 

with 

C ab , cd = r be f cde + f ade f cbe_ (1Q) 



In order to understand the importance of a non-constant background 
field, it is instructive to consider first the case of a homogeneous field A l Ax) = 
A\, with A\ ~ 0(l/g). For the case of time and space translation invariant 
correlators ([5]) we consider their Fourier transform F(p) and p(j>)o Similar 
to ([5]) we introduce statistical, 11(f), an d spectral, II(p), components of self- 
energies defined as [9] 



(P) 



L1 (F)ab 



nT,..(p) 



(p)a& V 









'(A)ab(P"> ' 



(11) 



where summation over repeated Lorentz and color indices is implied. The 
translation invariant propagators ([5]) and self-energies (jlip obey the iden- 
tity [15] 

nf; ) >)P 6 4(p)-nf; ) >)^>) = o, (12) 

which can be directly verified by plugging in the above definitions. Equa- 
tion (|12p is well-known in nonequilibrium physics and will be the starting 
point for our calculation. In the language of Boltzmann dynamics it states 
that "gain terms" equal "loss terms" for which stationarity is achieved [9]. 
Thermal equilibrium trivially solves (|12p . which we do not consider in the 
following. Instead, we will look for possible non-thermal scaling solutions. 

Decomposing the one-loop self-energy shown in Fig. 0] into its statistical 
(real) and spectral (imaginary) part, one obtains 

K )ef (p-,A) = £ f (27r)U^(p + q + k)V e ^(p, q ,k;A) 

^ J ok 



F^ a {q)F^{k) + - P h p a a {q)~p%{k) 



Vftdi-P'-Q'-k'-A) 



K)eM^) 



(27r)U^(p + q + k)V e ^(p,q,k;A) 



qk 



^dc I 



~ba 



pdc 



F?M~P%(k)+i% a {q)F%{k) V^(-p,-q,-k;A) (13) 



rv/36, 
'fbd 



with the notation f = f d 4 q/(27r) 4: and the symmetry property of the anti- 



commutator and commutator functions 



Kt(-p) 



F%b) 



P$A 



-p) 



-p%(p)- 



(14) 



x We introduce a —i in Fourier transforms of the spectral (p-) and retarded/advanced 
components, such as p(p) — — i f d 4 'xe lp ' J ' x p(x), while F(p) — fd 4 xe w > J ' x F(x). 
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Following similar steps as in Refs. |12|, 121 J . scaling solutions can be effi- 
ciently identified by integrating (|12p over external spatial momentum p and 
suitable scaling transformations for coordinates. In this way the problem 
can be reduced to simple algebraic conditions for scaling exponents. Non- 
thermal scaling solutions may be obtained in the classical regime, where 
expectation values of anti-commutators are much larger than commutators, 
F 2 S> p 2 ■ This is analogous to what is done in the context of weak Kol- 
mogorov wave turbulence using kinetic equations in the regime of sufficiently 
large occupation numbers [25]. In contrast, for lower occupancies of order 
one, dissipative or quantum corrections will obstruct scaling. 

The stationarity condition then reads 

= ^ / (2n)U( 4 \p + q + k)V^(p,q,k;A)V^ d 5 (-p,-q,-h,A) 

" Jpqk 

F b ^q)F s d c(k)pll(p) + F^q)^ 

> „ ' „ ' „ '. 

(I) (II) (III) (15) 

where J = J* d 3 p/(27r) 3 . We are looking for scaling solutions which behave 
as 

Kl(sp) = \s\-^F$(p) , pf v (sp) = \s\- 2 sgn(s) p* (p) (16) 

under rescaling with the real parameter s. This behavior reflects the scaling 
of the spectral function with the canonical dimension and takes into account 
a possible occupation number exponent k for the statistical function. This is 
equivalent to what is obtained from using a kinetic approach. We emphasize 
that no explicit gauge fixing has to be applied here and the calculation goes 
through for all gauges which admit scaling solutions of the form (|16|) . The 
scaling properties of the three-gluon vertex ((JJ) are different for the standard 
tree- level part (|Bj) and the background- field part (|SJ). To parametrize this 
behavior in a compact way, we write 

V^(sp, sq, sk) = s v V$(p, q, k) , (17) 

where v = 1 for the standard tree-level vertex and v = for the background 
field part. 

We now map the terms (II) and (III) in f|15|) using scaling transforma- 
tions such that they have the same form as (I) up to momentum dependent 
prefactors. First, we apply for (II) the transformation 

p p p 

q -)■ -^q , k -)■ -^p , p -► -jpk. (18) 

11 



The absolute value of the Jacobian for the frequency part of this transfor- 
mation is |p°/A; | 3 . The same procedure applies to the term (III) where the 
roles of (k, 7, 5, c, d) and (q, a, )3, a, b) are interchanged. Taking also into ac- 
count the symmetry of the three-gluon vertex under exchange of combined 
momentum, color and Lorentz index leads us to 







El 
2 



/ (2^) 4 6® (p + q + k) V e ^(p, q, k; A)V$(-p, -q, -k; A) 



Ipqk 

xFfc{q)F%(k)p&<p) 



1 + 



p u 



sgn 



k a 



+ 



p° 



Here 



A - 3 • 3 + 3^-4^-2(2 + k) v -2_W 1} W 2) 

VV 



• (19) 



(20) 



measure 5's FF 



P 



and v^- 1 ' , y( 2 > are the scaling exponents of the respective parts of the two 
vertex functions appearing in the one-loop self-energy 

As is well known, the above analysis in the presence of a homogeneous 
background field is complicated by the fact that to order g 2 the individual 
terms in (|15j) vanish, since the corresponding processes are kinematically 
forbidden on-shell, i.e. for p° = ±[p|. Only if the background field has 
a residual (space-) time dependence, A(t), the phase space is opened up 
by the associated frequency contribution of the field and it becomes kine- 
matically allowed. For the above scaling analysis, this amounts to taking 
iA 1 ) = t;( 2 ) = such that all contributions not involving the background field 
vanish. After taking this into account, the analysis is analogous to the dis- 
cussion in scalar field theory where the corresponding dynamics has been 
analyzed in great detail [161 02] ■ 

Therefore, it is also not surprising that the gauge theory in the presence 
of a time-varying background may exhibit the same scaling exponents as a 
scalar field theory with quartic self-interaction in the presence of a time- 
dependent condensate [Hj. Indeed, if we set A = —1 in (fT9|) then the 
5(p +q° + k°) in the integrand of (|19[) ensures the vanishing of (1 4- k° /p° + 
q°/p°), using \k°/p°\ sgn(p°/k°) = k°/p° with nonzero p°. From this one 
can directly read off with fj20f) the scaling solution 

A = -1 =* k = - (21) 

2 v ; 

associated to an energy cascade towards higher momenta [251 [T6| [T2| I2T] . 
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4 Conclusion 

In this work we have studied the nonequilibrium time evolution of SU(2) 
gauge theory starting from initial over-population. The classical-statistical 
simulations reveal a quick evolution towards an approximate scaling be- 
havior, with a subsequent quasi-stationary evolution well described by the 
scaling exponent k = 3/2 for a time duration of hundreds in units of Q s . 
These lattice results are compared to resummed perturbative estimates at 
not too low momenta, where the occupancies allow for a kinetic description. 
Remarkably, the value for the scaling exponent may be understood as aris- 
ing from the leading 0(g 2 ) contribution in the presence of a time-dependent 
background field. The phenomenon is associated to weak wave turbulence 
describing an energy cascade towards higher momenta. This particular as- 
pect is very similar to what is observed for scalar theories, where an effective 
cubic interaction arises because of the presence of a time-dependent Bose 
condensate. 

Of course, there are important differences between gauge theories and 
scalar theories. Competing elastic and inelastic scattering processes are 
expected to prevent a parametrically long time scale for the dominance 
of a background field dependent contribution in gauge theories, which is 
consistent with our findings. In particular, the interpretation of the most 
infrared modes in terms of distribution functions or even the notion of a 
Bose condensate in gauge theories is non-trivial. The question of whether 
an intermediate-time behavior is due to a condensate or to very soft modes 
associated to some slowly varying background-field in time and space may 
not be too critical after all in view of suggested potential signatures [26] , 
For this one also has to consider these questions in classical simulations 
including the relevant physics of longitudinal expansion in the context of 
heavy- ion collisions for sufficiently high Q s . 
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